Banded Dotterel Moult Study

Exploration of dataset

Authors

Affiliations

Luke Eberhart-Hertel

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Bashar Jarayseh

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Ailsa Howard

South Bay Banded Dotterel Project, Kaikoura, New Zealand

Tony Habraken

Port Waikato Banded Dotterel Project, Port Waikato, New Zealand

Emma Williams

Department of Conservation, Christchurch, New Zealand

Colin O`Donnell

Department of Conservation, Christchurch, New Zealand

Bart Kempenaers

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Published

January 22, 2025

Code

knitr::opts_chunk$set(cache = TRUE)

Prerequisites

R packages

Plotting themes

Visually explore the moult data extracted from photos

Note: data last modified 7-Nov-2024 This dataset contains the scored moult data for all usable photos from the 2021-2022, 2022-2023, and 2023-2024 breeding seasons at Kaikoura

Plot the sampling distribution of Ailsa’s photos across the seasons

2021-2022 season data

Code

# 2021-2022 season data
ggplot(data = dat %>% filter(season == 2021 & sex == "F") %>% mutate(score = as.numeric(score), 
                                                                     rings_comb_ = reorder(rings_comb, n_photos, decreasing = TRUE))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ rings_comb_, ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +
  xlab("week") +
  ggtitle("Females from the 2021-2022 breeding season in Kaikoura")

Code

ggplot(data = dat %>% filter(season == 2021 & sex == "M") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) +
  xlab("week") +
  ggtitle("Males from the 2021-2022 breeding season in Kaikoura")

2022-2023 season data

Code

# 2022-2023 season data (not as good coverage as the previous season)
ggplot(data = dat %>% filter(season == 2022 & sex == "F") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +
  xlab("week") +
  ggtitle("Females from the 2022-2023 breeding season in Kaikoura")

Code

ggplot(data = dat %>% filter(season == 2022 & sex == "M") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) +
  xlab("week") +
  ggtitle("Males from the 2022-2023 breeding season in Kaikoura")

2023-2024 season data

Code

# 2023-2024 season data 
ggplot(data = dat %>% filter(season == 2023 & sex == "F") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2023-07-01"), as.Date("2024-05-01"))) +
  xlab("week") +
  ggtitle("Females from the 2023-2024 breeding season in Kaikoura")

Code

ggplot(data = dat %>% filter(season == 2023 & sex == "M") %>% mutate(score = as.numeric(score))) +
  geom_point(aes(y = 1, x = date, fill = score), 
             pch = 21, color = "black", size = 3) +
  facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") +
  scale_fill_gradient(high = "#cc4c02", low = "white") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        strip.text.y.right = element_text(angle = 0)) +
  scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), 
               date_breaks = "3 week", 
               limits = c(as.Date("2023-07-05"), as.Date("2024-05-01"))) +
  xlab("week") +
  ggtitle("Males from the 2023-2024 breeding season in Kaikoura")

Statistical analysis of moult dynamics

Sex-specific differences in breeding plumage

sex-specific breeding plumage distributions

Code

ind_breeding_scores %>% 
  group_by(rings_comb, sex) %>% 
  summarise(max_breeding_score = max(max_breeding_score)) %>% 
  ungroup() %>% 
  group_by(sex, max_breeding_score) %>% 
  summarise(n = n()) %>% 
  bind_rows(., data.frame(sex = c("F", "M"), max_breeding_score = c(7, 3), n = c(0, 0))) %>%
  ungroup() %>% 
  ggplot() +
  geom_bar(aes(x = max_breeding_score, y = n, fill = sex), 
           alpha = 0.5, stat = "identity", position = "dodge", width = 0.5) +
  luke_theme +
  theme(legend.position = c(0.15, 0.9)) +
  xlab("Maximum individual rufous band score during breeding season") +
  ylab("Frequency") +
  scale_colour_brewer(palette = "Dark2", direction = -1,
                      name = "Sex",
                      labels = c("Female (N = 49)", "Male (N = 35)")) +
  scale_fill_brewer(palette = "Dark2", direction = -1,
                    name = "Sex",
                    labels = c("Female (N = 49)", "Male (N = 35)"))

effect-size table for sex-specific breeding plumage

Males tend to have a more full and immaculate rufous bands than females

Code

effects_mod_max_score <- 
  as.data.frame(allEffects(mod_max_score)$sex) %>%
  mutate(x_mean = ifelse(sex == "F", as.numeric(factor(sex)) + 0.2,
                         ifelse(sex == "M", as.numeric(factor(sex)) - 0.2, NA)),
         x_data = ifelse(sex == "F", as.numeric(factor(sex)), 
                          ifelse(sex == "M", as.numeric(factor(sex)), NA)))


ggplot() +
  geom_errorbar(data = effects_mod_max_score, 
                aes(x = x_mean, ymin = lower, ymax = upper), width = 0.05) +
  geom_point(data = effects_mod_max_score, 
             aes(x = x_mean, y = fit, fill = sex), shape = 21, size = 5) +
  geom_jitter(data = ind_breeding_scores,
             aes(x = as.numeric(factor(sex)), 
                 y = max_breeding_score, fill = sex), 
             width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75) +
  theme(legend.position = "none",
        text = element_text(family = "lato"),
        axis.ticks = element_blank(),
        axis.text = element_text(size = 10, colour = "grey40"),
        axis.title.y = element_text(size = 11, colour = "grey40"),
        axis.text.y = element_text(size = 10, colour = "grey40"),
        plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
        panel.spacing = unit(1, "cm"),
        panel.background = element_rect(fill = 'transparent', color = NA),
        plot.background = element_rect(fill = 'transparent', color = NA), 
        strip.background = element_rect(fill = 'grey90', color = NA),
        strip.text = element_text(size = 13, colour = "grey40"),
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        axis.title.x = element_blank()) +
  scale_x_continuous(limits = c(0.8, 2.2), 
                     breaks = c(1.1, 1.9), 
                     labels = c("female", "male")) +
  ylab("maximum breeding plumage score observed") +
  scale_fill_brewer(palette = "Set2") +
  scale_y_continuous(limits = c(3.8, 7.2))

Age- and sex-specific variation in breeding plumage (of known-aged birds)

effect-size table for sex- and age-specific variation in breeding plumage

no effect of age on the maximum extent of an individual’s rufous band score while controlling for underlying sex-specific variation

Migratory status- and sex-specific variation in breeding plumage (of known-status birds)

Code

# draw gt table
mod_max_score_table_status_sex <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_label(comp_name = html("<i>Banded Dotterel rufous band score</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = c(estimate),
             rows = 1:10,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = c(estimate),
             rows = 11:13,
             decimals = 0,
             use_seps = FALSE) %>% 
  sub_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = c(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(80))

effect-size table and plot for sex- and migratory-status breeding plumage

no effect of migratory status on an individual’s rufous band score while controlling for underlying sex-specific variation

Within-individual moult dynamics across sexes

Code

# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in sex while accounting for sex differences in maximum rufous band score.
ind_prop_molt_scores <- 
  dat %>% 
  mutate(score = as.numeric(score)) %>% 
  left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% 
  filter(!is.na(max_breeding_score)) %>% 
  mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
  select(rings_comb, season, date_J, sex, score, max_breeding_score, prop_molt_score) %>% 
  arrange(rings_comb, season, date_J) %>% 
  distinct() %>% 
  filter(!is.na(prop_molt_score)) %>%
  group_by(sex) %>%
  mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>%
  ungroup()

# Assess sample sizes of each sex
ind_prop_molt_scores %>% 
  group_by(sex) %>% 
  summarise(n_distinct(rings_comb))

# A tibble: 2 × 2
  sex   `n_distinct(rings_comb)`
  <chr>                    <int>
1 F                           49
2 M                           35

Code

# mixed effects binomial model comparing sex and date effect on the changes in moult scores while controlling for between sex differences in the maximum breeding score
mod1 <- 
  lme4::glmer(prop_molt_score ~ 
                date_J * sex + std_max_breeding_score +
         (1 | rings_comb),
       data = ind_prop_molt_scores, 
       family = binomial,
       control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

# strong date effect, but no sex effect
tbl_regression(mod1, intercept = TRUE, 
               label = list(date_J ~ "Date", sex ~ "Sex"))

Characteristic	log(OR)¹	95% CI¹	p-value
(Intercept)	20	15, 24	<0.001
Date	-0.09	-0.11, -0.07	<0.001
Sex
F	—	—
M	1.1	-5.7, 8.0	0.8
std_max_breeding_score	-0.29	-0.65, 0.07	0.11
Date * Sex
Date * M	-0.01	-0.04, 0.03	0.7
¹ OR = Odds Ratio, CI = Confidence Interval

Code

# extract predicted trends
pr <- ggeffects::predict_response(mod1, c("date_J [30:293]", "sex"))

# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

# join the back-transformed dates to model fits
mod1_fits <- 
  as.data.frame(pr) %>% 
  rename(date_J = x,
         sex = group) %>% 
  left_join(., dates_for_plot, by = "date_J")

ind_prop_molt_scores <- 
  ind_prop_molt_scores %>% 
  left_join(., dates_for_plot, by = "date_J")

# plot the model
ggplot() +
  geom_line(data = mod1_fits, 
            aes(x = date, y = predicted, color = sex)) +
  geom_ribbon(data = mod1_fits, 
              aes(x = date, ymax = conf.high, ymin = conf.low, fill = sex),
              lwd = 1, alpha = 0.25) +
  geom_jitter(data = ind_prop_molt_scores,
              aes(x = date, y = prop_molt_score, color = sex), 
              alpha = 0.25, shape = 20, width = 0, height = 0.025) +
  theme(legend.position = c(0.3, 0.2),
        legend.justification = c(1, 0),
        text = element_text(family = "lato"),
        axis.ticks = element_blank(),
        axis.text.x = element_text(size = 10, 
                                   angle = 45, 
                                   hjust = 1, 
                                   vjust = 1),
        axis.text.y = element_text(size = 10, colour = "grey40"),
        axis.title.y = element_text(size = 11, colour = "grey40"),
        plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
        panel.spacing = unit(1, "cm"),
        panel.background = element_rect(fill = 'transparent', color = NA),
        plot.background = element_rect(fill = 'transparent', color = NA), 
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        axis.title.x = element_blank()) +
  ylab("proportion intact of individual-based maximum rufous band") +
  scale_x_date(date_labels = "%B", 
               expand = c(0.01, 0.01), 
               date_breaks = "1 month") +
  scale_colour_brewer(palette = "Dark2", direction = -1,
                      name = "Sex",
                      labels = c("Female (N = 49)", "Male (N = 35)")) +
  scale_fill_brewer(palette = "Dark2", direction = -1,
                    name = "Sex",
                    labels = c("Female (N = 49)", "Male (N = 35)"))

Code

# wrangle the dates of plumage retention and moult for subsequent analysis
# Step 1: Calculate the max_score for each individual and season
max_scores <- ind_prop_molt_scores %>%
  group_by(rings_comb, season) %>%
  summarise(max_score = max(prop_molt_score, na.rm = TRUE), .groups = "drop")

# Step 2: Merge the max_scores back into the original data
mid_point_dates <- 
  ind_prop_molt_scores %>%
  left_join(max_scores, by = c("rings_comb", "season")) %>%  # Add the max_score
  group_by(rings_comb, season) %>%
  arrange(date_J) %>%  # Ensure data is sorted by date_J

  # Step 3: Find the latest date_J for the maximum prop_molt_score
  filter(prop_molt_score == max_score) %>%
  slice_max(order_by = date_J, n = 1, with_ties = FALSE) %>%
  rename(latest_max_date_J = date_J) %>%
  
  # Step 4: Find the earliest date_J where prop_molt_score < max_score
  right_join(
    ind_prop_molt_scores %>%
      left_join(max_scores, by = c("rings_comb", "season")) %>%  # Add the max_score here too
      group_by(rings_comb, season) %>%
      arrange(date_J) %>%
      filter(prop_molt_score < max_score) %>%
      slice_min(order_by = date_J, n = 1, with_ties = FALSE) %>%
      rename(earliest_decline_date_J = date_J),
    by = c("rings_comb", "season")
  ) %>%
  
  # Step 5: Compute the mid-point date_J
  mutate(mid_point_date_J = round((latest_max_date_J + earliest_decline_date_J) / 2)) %>%
  
  # Select relevant columns
  select(rings_comb, season, latest_max_date_J, earliest_decline_date_J, mid_point_date_J) %>%
  distinct()

dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

moult_commence_ring_season <- 
  left_join(mid_point_dates, 
            dates_for_plot %>% rename(mid_point_date_J = date_J), 
            by = "mid_point_date_J")

Within-individual moult dynamics across migratory status

Code

# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status.
ind_prop_molt_scores <- 
  dat %>% 
  mutate(score = as.numeric(score)) %>% 
  left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% 
  filter(!is.na(max_breeding_score)) %>% 
  mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
  select(rings_comb, season, date_J, sex, migratory_status, score, max_breeding_score, prop_molt_score) %>% 
  arrange(rings_comb, season, date_J) %>% 
  distinct() %>% 
  filter(!is.na(prop_molt_score)) %>% 
  filter(migratory_status %in% c("R", "M"))

# Assess sample sizes of each migratory status
ind_prop_molt_scores %>% 
  group_by(migratory_status) %>% 
  summarise(n_distinct(rings_comb))

# A tibble: 2 × 2
  migratory_status `n_distinct(rings_comb)`
  <chr>                               <int>
1 M                                       5
2 R                                      16

Code

# mixed effects binomial model comparing migratory status and date effect on the changes in moult scores
mod1_status <- 
  lme4::glmer(prop_molt_score ~ 
                date_J * migratory_status +
         (1 | rings_comb),
       data = ind_prop_molt_scores, 
       family = binomial,
       control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

# strong date effect, but no status effect
tbl_regression(mod1_status, intercept = TRUE, 
               label = list(date_J ~ "Date", migratory_status ~ "Migratory status"))

Characteristic	log(OR)¹	95% CI¹	p-value
(Intercept)	25	10, 40	0.001
Date	-0.11	-0.19, -0.04	0.002
Migratory status
M	—	—
R	12	-17, 41	0.4
Date * Migratory status
Date * R	-0.06	-0.20, 0.08	0.4
¹ OR = Odds Ratio, CI = Confidence Interval

Code

# extract predicted trends
pr <- ggeffects::predict_response(mod1_status, c("date_J [30:293]", "migratory_status"))

# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

# join the back-transformed dates to model fits
mod1_fits <- 
  as.data.frame(pr) %>% 
  rename(date_J = x,
         migratory_status = group) %>% 
  left_join(., dates_for_plot, by = "date_J")

ind_prop_molt_scores <- 
  ind_prop_molt_scores %>% 
  left_join(., dates_for_plot, by = "date_J")

# plot the model
ggplot() +
  geom_line(data = mod1_fits, 
            aes(x = date, y = predicted, color = migratory_status)) +
  geom_ribbon(data = mod1_fits, 
              aes(x = date, ymax = conf.high, ymin = conf.low, fill = migratory_status),
              lwd = 1, alpha = 0.25) +
  geom_jitter(data = ind_prop_molt_scores,
              aes(x = date, y = prop_molt_score, color = migratory_status), 
              alpha = 0.25, shape = 20, width = 0, height = 0.025) +
  theme(legend.position = c(0.3, 0.2),
        legend.justification = c(1, 0),
        text = element_text(family = "lato"),
        axis.ticks = element_blank(),
        axis.text.x = element_text(size = 10, 
                                   angle = 45, 
                                   hjust = 1, 
                                   vjust = 1),
        axis.text.y = element_text(size = 10, colour = "grey40"),
        axis.title.y = element_text(size = 11, colour = "grey40"),
        plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
        panel.spacing = unit(1, "cm"),
        panel.background = element_rect(fill = 'transparent', color = NA),
        plot.background = element_rect(fill = 'transparent', color = NA), 
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        axis.title.x = element_blank()) +
  ylab("proportion intact of individual-based maximum rufous band") +
  scale_x_date(date_labels = "%B", 
               expand = c(0.01, 0.01), 
               date_breaks = "1 month") +
  scale_colour_brewer(palette = "Set1", direction = -1,
                      name = "Migratory status",
                      labels = c("Resident (N = 16)", "Migratory (N = 5)")) +
  scale_fill_brewer(palette = "Set1", direction = -1,
                    name = "Migratory status",
                      labels = c("Resident (N = 16)", "Migratory (N = 5)"))

Within-individual moult dynamics across ages

Code

# calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation between age groups.
ind_prop_molt_scores <- 
  dat %>% 
  mutate(score = as.numeric(score)) %>% 
  left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% 
  filter(!is.na(max_breeding_score)) %>% 
  mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
  select(rings_comb, season, date_J, sex, age, score, max_breeding_score, prop_molt_score) %>% 
  arrange(rings_comb, season, date_J) %>% 
  distinct() %>% 
  filter(!is.na(prop_molt_score)) %>% 
  filter(!is.na(age))

# Assess sample sizes of each sex
ind_prop_molt_scores %>% 
  group_by(age) %>% 
  summarise(n_distinct(rings_comb))

# A tibble: 6 × 2
    age `n_distinct(rings_comb)`
  <dbl>                    <int>
1     1                        6
2     2                        5
3     3                       12
4     4                       13
5     5                        4
6     6                        2

Code

# mixed effects binomial model comparing sex and date effect on the changes in moult scores
mod1_age <- 
  lme4::glmer(prop_molt_score ~ 
                date_J + age +
         # date_J * sex +
         (1 | rings_comb),# + (1 | season),
       data = ind_prop_molt_scores, 
       family = binomial)#,
       # control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06)))

plot(allEffects(mod1_age))

Code

# strong date effect, but no status effect
tbl_regression(mod1_age, intercept = TRUE, 
               label = list(date_J ~ "Date", age ~ "Age"))

Characteristic	log(OR)¹	95% CI¹	p-value
(Intercept)	53	18, 88	0.003
Date	-0.28	-0.47, -0.10	0.003
Age	2.2	0.61, 3.7	0.006
¹ OR = Odds Ratio, CI = Confidence Interval

Code

# extract predicted trends
pr <- ggeffects::predict_response(mod1_age, c("date_J [30:293]", "age"))

# back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date
dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

# join the back-transformed dates to model fits
mod1_fits <- 
  as.data.frame(pr) %>% 
  rename(date_J = x,
         age = group) %>% 
  left_join(., dates_for_plot, by = "date_J")

ind_prop_molt_scores <- 
  ind_prop_molt_scores %>% 
  left_join(., dates_for_plot, by = "date_J")

# plot the model
ggplot() +
  geom_line(data = mod1_fits, 
            aes(x = date, y = predicted, color = age)) +
  geom_ribbon(data = mod1_fits, 
              aes(x = date, ymax = conf.high, ymin = conf.low, fill = age),
              lwd = 1, alpha = 0.25) +
  geom_jitter(data = ind_prop_molt_scores,
              aes(x = date, y = prop_molt_score, color = as.factor(age)),
              alpha = 0.25, shape = 20, width = 0, height = 0.025) +
  theme(legend.position = c(0.3, 0.2),
        legend.justification = c(1, 0),
        text = element_text(family = "lato"),
        axis.ticks = element_blank(),
        axis.text.x = element_text(size = 10, 
                                   angle = 45, 
                                   hjust = 1, 
                                   vjust = 1),
        axis.text.y = element_text(size = 10, colour = "grey40"),
        axis.title.y = element_text(size = 11, colour = "grey40"),
        plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
        panel.spacing = unit(1, "cm"),
        panel.background = element_rect(fill = 'transparent', color = NA),
        plot.background = element_rect(fill = 'transparent', color = NA), 
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        axis.title.x = element_blank()) +
  ylab("proportion intact of individual-based maximum rufous band") +
  scale_x_date(date_labels = "%B", 
               expand = c(0.01, 0.01), 
               date_breaks = "1 month", 
               limits = c(as.Date("2021-10-01"), as.Date("2022-05-01"))) +
  scale_colour_brewer(palette = "Spectral", direction = -1,
                      name = "Age",
                      labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", 
                                 "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) +
  scale_fill_brewer(palette = "Spectral", direction = -1,
                    name = "Age",
                      labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", 
                                 "4 (N = 13", "5 (N = 4)", "6 (N = 2)"))

Breeding status analysis

Code

# explore datasets provided by Ted
## 2021 season
# import and consolidate key columns
breeding_data_2021 <- 
  read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx", 
                     sheet = "Kaikoura_NestVisit2021_1", 
                     col_types = "text") %>%
  mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% 
  mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% 
  mutate(date = as.Date(visit_date_nz)) %>% 
  select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% 
  rename(nest_fate = `Nest Status`,
         eggs = `Egg count`,
         hatched = `Number Hatched`,
         fledged = `Number Fledged`,
         chicks = `Number of chicks seen`)

# Define a function to extract the desired patterns and ensure enough columns
extract_patterns <- function(text) {
  # Replace vertical bar '|' with an empty string
  text <- gsub("\\|", "", text)
  
  # Extract all 4 or 5-character long texts that start with R or r
  matches_r <- str_extract_all(text, "\\b[Rr]\\w{3,4}\\b")[[1]]
  
  # Extract UB, UN, UBF, UBM as standalone patterns
  matches_ub_un <- str_extract_all(text, "\\b[Uu][BbNnFfMm]\\b")[[1]]
  
  # Combine matches
  matches <- c(matches_r, matches_ub_un)
  
  # Ensure all matches are capitalized
  matches <- toupper(matches)
  
  # Ensure there are exactly two columns, filling with NA if necessary
  result <- c(matches, rep(NA, 2 - length(matches)))
  
  # Return the result as a named vector
  return(setNames(result, c("parent1", "parent2")))
}

# Apply the function to the text column and store the results in new columns
extracted_parents <- t(sapply(breeding_data_2021$NestID, extract_patterns))

breeding_data_2021 <-
  bind_cols(breeding_data_2021, extracted_parents) %>% 
  select(-parent2) %>% 
  rename(parent = parent1) %>% 
  bind_rows(bind_cols(breeding_data_2021, extracted_parents) %>% 
              select(-parent1) %>% 
              rename(parent = parent2)) %>% 
  filter(!is.na(parent)) %>% 
  filter(parent %in% (dat %>% filter(season == 2021) %>% pull(rings_comb) %>% unique())) %>%
  filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>%
  group_by(parent) %>%
  arrange(desc(date)) %>%
  slice(1) %>%
  ungroup() %>%
  arrange(desc(date)) %>%
  mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>%
  rename(last_date_breeding = date) %>% 
  select(parent, last_date_breeding, nest_id, OBJECTID)

# import and consolidate key columns
breeding_data_2021_ <-
  read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx", 
                     sheet = "Kaikoura_DotterelNest2021_0", 
                     col_types = "text") %>%
  mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"),
         fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"),
         found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% 
  mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"),
         fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"),
         found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% 
  mutate(hatch_date = as.Date(hatch_date_nz),
         fail_date = as.Date(fail_date_nz),
         found_date = as.Date(found_date_nz)) %>%
  mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% 
  select(OBJECTID, nest_id, found_date, hatch_date, fail_date)

breeding_2021_final <- left_join(breeding_data_2021, breeding_data_2021_, by = "nest_id") %>%
  mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | 
                               ((!is.na(hatch_date) & found_date < hatch_date)) | 
                               ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) |
                               ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% 
  group_by(parent) %>% 
  mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% 
  mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% 
  arrange(desc(days_diff)) %>% 
  mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>% 
  arrange(desc(last_date_breeding_final)) %>% 
  select(parent, last_date_breeding_final) %>% 
  # specify the season as the first calender year
  mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final)))

###########

## 2022 season

# import and consolidate key columns
breeding_data_2022 <- 
  read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx", 
                     sheet = "Kaikoura_NestVisit_1", 
                     col_types = "text") %>%
  mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% 
  mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% 
  mutate(date = as.Date(visit_date_nz)) %>% 
  select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% 
  rename(nest_fate = `Nest Status`,
         eggs = `Egg count`,
         hatched = `Number Hatched`,
         fledged = `Number Fledged`,
         chicks = `Number of chicks seen`)


# Apply the function to the text column and store the results in new columns
extracted_parents2 <- t(sapply(breeding_data_2022$NestID, extract_patterns))

breeding_data_2022 <-
  bind_cols(breeding_data_2022, extracted_parents2) %>% 
  select(-parent2) %>% 
  rename(parent = parent1) %>% 
  bind_rows(bind_cols(breeding_data_2022, extracted_parents2) %>% 
              select(-parent1) %>% 
              rename(parent = parent2)) %>% 
  filter(!is.na(parent)) %>% 
  filter(parent %in% (dat %>% filter(season == 2022) %>% pull(rings_comb) %>% unique())) %>%
  filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>%
  group_by(parent) %>%
  arrange(desc(date)) %>%
  slice(1) %>%
  ungroup() %>%
  arrange(desc(date)) %>%
  mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>%
  rename(last_date_breeding = date) %>% 
  select(parent, last_date_breeding, nest_id, OBJECTID)

# import and consolidate key columns
breeding_data_2022_ <-
  read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx", 
                     sheet = "Kaikoura_DotterelNest_0", 
                     col_types = "text") %>%
  mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"),
         fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"),
         found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% 
  mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"),
         fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"),
         found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% 
  mutate(hatch_date = as.Date(hatch_date_nz),
         fail_date = as.Date(fail_date_nz),
         found_date = as.Date(found_date_nz)) %>%
  mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% 
  select(OBJECTID, nest_id, found_date, hatch_date, fail_date)

breeding_2022_final <- left_join(breeding_data_2022, breeding_data_2022_, by = "nest_id") %>%
  mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | 
                               ((!is.na(hatch_date) & found_date < hatch_date)) | 
                               ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) |
                               ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% 
  group_by(parent) %>% 
  mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% 
  mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% 
  arrange(desc(days_diff)) %>% 
  mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>% 
  arrange(desc(last_date_breeding_final)) %>% 
  select(parent, last_date_breeding_final) %>% 
  # specify the season as the first calender year
  mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final)))



###########

## 2023 season

# import and consolidate key columns
breeding_data_2023 <- 
  read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx", 
                     sheet = "Kaikoura_NestVisit_1", 
                     col_types = "text") %>%
  mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% 
  mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% 
  mutate(date = as.Date(visit_date_nz)) %>% 
  select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% 
  rename(nest_fate = `Nest Status`,
         eggs = `Egg count`,
         hatched = `Number Hatched`,
         fledged = `Number Fledged`,
         chicks = `Number of chicks seen`)


# Apply the function to the text column and store the results in new columns
extracted_parents3 <- t(sapply(breeding_data_2023$NestID, extract_patterns))

breeding_data_2023 <-
  bind_cols(breeding_data_2023, extracted_parents3) %>% 
  select(-parent2) %>% 
  rename(parent = parent1) %>% 
  bind_rows(bind_cols(breeding_data_2023, extracted_parents3) %>% 
              select(-parent1) %>% 
              rename(parent = parent2)) %>% 
  filter(!is.na(parent)) %>% 
  filter(parent %in% (dat %>% filter(season == 2023) %>% pull(rings_comb) %>% unique())) %>%
  filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>%
  group_by(parent) %>%
  arrange(desc(date)) %>%
  slice(1) %>%
  ungroup() %>%
  arrange(desc(date)) %>%
  mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>%
  rename(last_date_breeding = date) %>% 
  select(parent, last_date_breeding, nest_id, OBJECTID)

# import and consolidate key columns
breeding_data_2023_ <-
  read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx", 
                     sheet = "Kaikoura_DotterelNest_0", 
                     col_types = "text") %>%
  mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"),
         fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"),
         found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% 
  mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"),
         fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"),
         found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% 
  mutate(hatch_date = as.Date(hatch_date_nz),
         fail_date = as.Date(fail_date_nz),
         found_date = as.Date(found_date_nz)) %>%
  mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% 
  select(OBJECTID, nest_id, found_date, hatch_date, fail_date)

breeding_2023_final <- left_join(breeding_data_2023, breeding_data_2023_, by = "nest_id") %>%
  mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | 
                               ((!is.na(hatch_date) & found_date < hatch_date)) | 
                               ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) |
                               ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% 
  group_by(parent) %>% 
  mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% 
  mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% 
  arrange(desc(days_diff)) %>% 
  mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>% 
  arrange(desc(last_date_breeding_final)) %>% 
  select(parent, last_date_breeding_final) %>% 
  # specify the season as the first calender year
  mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final)))


# combine breeding data from the 3 seasons into one table

breeding_data_all <- bind_rows(breeding_2021_final,breeding_2022_final,breeding_2023_final) %>% 
  rename(rings_comb = parent) %>% 
  # change to Julian date shifted for the Southern Hemisphere (1 = July 1)
  mutate(breeding_date_J = as.numeric(format(last_date_breeding_final + 181, "%j"))) %>% 
  ungroup() %>% 
  # standardize the last breeding date by year
  mutate(breeding_date_J_s = scale_by(breeding_date_J ~ season, ., scale = 0)) %>% 
  # make the scaled date variable numeric class
  mutate(breeding_date_J_snum = as.numeric(breeding_date_J_s))

Relationships between breeding schedule and plumage

Code

ind_prop_molt_scores_breeding_date <-
  dat %>% 
  mutate(score = as.numeric(score)) %>% 
  left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% 
  filter(!is.na(max_breeding_score)) %>% 
  mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>%
  select(rings_comb, season, date_J, sex, migratory_status, age, score, max_breeding_score, prop_molt_score) %>% 
  arrange(rings_comb, season, date_J) %>% 
  distinct() %>% 
  filter(!is.na(prop_molt_score)) %>% 
  full_join(breeding_data_all, by = c("rings_comb", "season")) %>% 
  mutate(date_diff = abs(breeding_date_J - date_J)) %>% 
  arrange(rings_comb, season, date_diff) %>% 
  group_by(rings_comb, season) %>% slice(1) %>% filter(date_diff < 30) %>% 
  dplyr::select(rings_comb, season, sex, breeding_date_J, breeding_date_J_snum, prop_molt_score, max_breeding_score, age) %>% 
  distinct() %>%
  group_by(sex) %>%
  mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>%
  ungroup() %>% 
  mutate(breeding_date_J_s = scale(breeding_date_J)) %>% 
  left_join(., moult_commence_ring_season, by = c("rings_comb", "season"))

effect of breeding plumage score on prolonged breeding

do individuals with a fuller breeding plumage (for their sex) breed later than the rest of the population? No effect

Code

mod_max_score_breeding_date <- 
  lme4::lmer(breeding_date_J_snum ~ 
                 std_max_breeding_score +
         (1 | rings_comb),
       data = ind_prop_molt_scores_breeding_date)

tbl_regression(mod_max_score_breeding_date, intercept = TRUE)

Characteristic	Beta	95% CI¹
(Intercept)	4.5	-0.64, 9.7
std_max_breeding_score	2.5	-2.7, 7.7
¹ CI = Confidence Interval

Code

plot(allEffects(mod_max_score_breeding_date), type = "response")

effect of breeding plumage score on prolonged breeding

do individuals with a late breeding attempts tend to have a fuller breeding plumage? No effect

Code

mod_breeding_date_max_score <- 
  lme4::lmer(std_max_breeding_score ~ 
                 breeding_date_J_snum +
         (1 | rings_comb),
       data = ind_prop_molt_scores_breeding_date)

summary(mod_max_score_breeding_date)

Linear mixed model fit by REML ['lmerMod']
Formula: breeding_date_J_snum ~ std_max_breeding_score + (1 | rings_comb)
   Data: ind_prop_molt_scores_breeding_date

REML criterion at convergence: 1075.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7013 -0.6280  0.1192  0.6279  1.9243 

Random effects:
 Groups     Name        Variance Std.Dev.
 rings_comb (Intercept)   0.0     0.0    
 Residual               795.4    28.2    
Number of obs: 114, groups:  rings_comb, 63

Fixed effects:
                       Estimate Std. Error t value
(Intercept)               4.539      2.641   1.718
std_max_breeding_score    2.475      2.665   0.929

Correlation of Fixed Effects:
            (Intr)
std_mx_brd_ 0.000 
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Code

plot(allEffects(mod_breeding_date_max_score), type = "response")

Relationships between breeding schedule and moult

effect of prolonged breeding on moult

do individuals with late breeding attempts have a delayed moult? Yes: individuals that breed later than the rest of the population in a given year tend to commence their complete/post-breeding moult later

Code

mod_breeding_date_prop_molt_mid <- 
  lme4::lmer(mid_point_date_J ~ 
                 breeding_date_J_snum +
         (1 | rings_comb),
       data = ind_prop_molt_scores_breeding_date)

tbl_regression(mod_breeding_date_prop_molt_mid, intercept = TRUE)

Characteristic	Beta	95% CI¹
(Intercept)	153	145, 160
breeding_date_J_snum	0.31	0.03, 0.58
¹ CI = Confidence Interval

Code

plot(allEffects(mod_breeding_date_prop_molt_mid), type = "response")

effect of prolonged breeding on prolonged retention of full plumage

do individuals with late breeding attempts retain their full breeding plumage longer? Yes: individuals that breed later than the rest of the population in a given year tend to retain their full breeding plumage longer

Code

mod_breeding_date_prop_molt_latest <- 
  lm(latest_max_date_J ~ breeding_date_J_snum,
       data = ind_prop_molt_scores_breeding_date)

tbl_regression(mod_breeding_date_prop_molt_latest, intercept = TRUE)

Characteristic	Beta	95% CI¹	p-value
(Intercept)	148	141, 155	<0.001
breeding_date_J_snum	0.59	0.34, 0.85	<0.001
¹ CI = Confidence Interval

Code

plot(allEffects(mod_breeding_date_prop_molt_latest), type = "response")

Relationships between age and extended breeding period

Code

# wrangle data to get ages and dates of last breeding for each season of all 
# birds available
age_breeding_date <-
  dat %>% 
  select(rings_comb, season, age, sex) %>% 
  full_join(breeding_data_all, by = c("rings_comb", "season")) %>% 
  distinct()

effect of age on prolonged breeding season

do older individuals conclude their breeding season later? No: there is strong evidence that younger individuals tend to have the latest breeding attempts in the population for a given year.

Code

mod_age_breeding_date <- 
  lme4::lmer(breeding_date_J_snum ~ age +
         (1 | rings_comb),
       data = age_breeding_date)

summary(mod_age_breeding_date)

Linear mixed model fit by REML ['lmerMod']
Formula: breeding_date_J_snum ~ age + (1 | rings_comb)
   Data: age_breeding_date

REML criterion at convergence: 312.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.3355 -0.5314  0.1108  0.5764  2.2836 

Random effects:
 Groups     Name        Variance Std.Dev.
 rings_comb (Intercept)    0      0.0    
 Residual               1069     32.7    
Number of obs: 33, groups:  rings_comb, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)   25.976     14.414   1.802
age           -7.608      3.468  -2.194

Correlation of Fixed Effects:
    (Intr)
age -0.919
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')

Code

plot(allEffects(mod_age_breeding_date), type = "response")

Code

# strong date effect, but no status effect
tbl_regression(mod_age_breeding_date, intercept = TRUE, 
               label = list(age ~ "Age"))

Characteristic	Beta	95% CI¹
(Intercept)	26	-2.3, 54
Age	-7.6	-14, -0.81
¹ CI = Confidence Interval

Code

# extract predicted trends
# I tried the whole range for the breeding timing (breeding_date_J [72:221]), but then decided to take some values with 20 days interval to visualize it better
pr <- ggeffects::predict_response(mod_age_breeding_date, c("age"))

dates_for_plot <-
  data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             breeding_date_ = c(1:365))
CL_dates_for_plot <-
  data.frame(conf.low_date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             conf.low_ = c(1:365))
CH_dates_for_plot <-
  data.frame(conf.high_date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             conf.high_ = c(1:365))

# join the back-transformed dates to model fits
mod1_fits <- 
  as.data.frame(pr) %>% 
  rename(age = x,
         # sex = group,
         breeding_date = predicted) %>% 
  mutate(breeding_date_ = round(breeding_date),
         conf.low_ = round(conf.low),
         conf.high_ = round(conf.high)) %>% 
  left_join(., dates_for_plot, by = "breeding_date_") %>% 
  left_join(., CL_dates_for_plot, by = "conf.low_") %>% 
  left_join(., CH_dates_for_plot, by = "conf.high_")
  
# plot the model
ggplot() +
  geom_line(data = mod1_fits, 
            aes(x = age, y = date)) +
  geom_ribbon(data = mod1_fits, 
              aes(x = age, ymax = conf.high_date, ymin = conf.low_date),
              lwd = 1, alpha = 0.25) +
  luke_theme +
  # theme(legend.position = c(0.3, 0.2),
  #       legend.justification = c(1, 0),
  #       strip.background = element_blank(),
  #       axis.title.x = element_blank(),
  #       axis.text.x = element_text(size = 10, 
  #                                  angle = 45, 
  #                                  hjust = 1, 
  #                                  vjust = 1)) +
  ylab("date of last breeding")

Meeting notes:

13-Nov-24: Bart, Luke, and Bashar via zoom

add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted.
drop pre-October data from the molt timing analysis.
next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.

--- title: "Banded Dotterel Moult Study" subtitle: | Exploration of dataset date: "`r format(Sys.time(), '%d %B, %Y')`" author: - name: Luke Eberhart-Hertel orcid: 0000-0001-7311-6088 email: luke.eberhart@bi.mpg.de url: https://www.bi.mpg.de/person/115852/2867 affiliations: - ref: bk - name: Bashar Jarayseh affiliations: - ref: bk - name: Ailsa Howard affiliations: - ref: ah - name: Tony Habraken affiliations: - ref: th - name: Emma Williams affiliations: - ref: ew - name: Colin O`Donnell affiliations: - ref: ew - name: Bart Kempenaers affiliations: - ref: bk affiliations: - id: bk number: 1 name: Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany - id: ah number: 2 name: South Bay Banded Dotterel Project, Kaikoura, New Zealand - id: th number: 2 name: Port Waikato Banded Dotterel Project, Port Waikato, New Zealand - id: ew number: 3 name: Department of Conservation, Christchurch, New Zealand format: html: toc: true code-fold: true code-tools: true self-contained: true highlight-style: github theme: Cosmo execute: warning: false cache: true editor_options: chunk_output_type: console --- ```{r} knitr::opts_chunk$set(cache = TRUE) ``` ## Prerequisites ### R packages ```{r, message=FALSE, results='hide', warning=FALSE, results='hide', cache=FALSE, include=FALSE} # - The following packages are needed for analysis and can be easily installed from [CRAN](http://cran.r-project.org/) or GitHub by running the following code chunk: # a vector of all the packages needed in the project packages_required_in_project <- c("tidyverse", "readxl", "RMark", "RColorBrewer", "patchwork", "mapview", "lubridate", "extrafont", "here", "DT", "leaflet", "sf", "leafpop", "tsibble", "corrplot", "gghalves", "gam", "pscl", "gamlss", "gt", "lme4", "ggpattern", "gtsummary", "effects", "lattice", "rptR", "partR2", "broom.mixed", "forcats", "glmmTMB", "sysfonts", "showtext", "patchwork", "standardize") # of the required packages, check if some need to be installed new.packages <- packages_required_in_project[!(packages_required_in_project %in% installed.packages()[,"Package"])] # install all packages that are not locally available if(length(new.packages)) install.packages(new.packages) # load all the packages into the current R session lapply(packages_required_in_project, require, character.only = TRUE) # set the home directory to where the project is locally based (i.e., to find # the relevant datasets to import, etc. here::set_here() ``` ### Plotting themes ```{r, message=FALSE, results='hide', warning=FALSE, include=FALSE} # - The following plotting themes, colors, and typefaces are used throughout the project: # Find fonts from computer that you want. Use regular expressions to do this # For example, load all fonts that are 'verdana' or 'Verdana' extrafont::font_import(pattern = "[V/v]erdana", prompt = FALSE) # check which fonts were loaded extrafont::fonts() extrafont::fonttable() extrafont::loadfonts() # load these into R # Add Google font Lato font_add_google("Lato", "lato") # Automatically use showtext to render text showtext_auto() # define the plotting theme to be used in subsequent ggplots luke_theme <- theme_bw() + theme( text = element_text(family = "Verdana"), legend.title = element_text(size = 10), legend.text = element_text(size = 8), axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 8), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 8), strip.text = element_text(size = 10), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_line(linewidth = 0.5, colour = "grey40"), axis.ticks.length = unit(0.2, "cm"), panel.border = element_rect(linetype = "solid", colour = "grey"), legend.position.inside = c(0.1, 0.9) ) ``` ## Visually explore the moult data extracted from photos Note: data last modified 7-Nov-2024 This dataset contains the scored moult data for all usable photos from the 2021-2022, 2022-2023, and 2023-2024 breeding seasons at Kaikoura ```{r, include = FALSE} # import data with all columns as character (so that no auto-formatting is done by R) dat <- read.csv("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Dataset_Molt_Final - modified 7.11.2024.csv", colClasses = "character") %>% mutate(Rings_comb = str_sub(Rings_comb, 1, 4)) %>% rename_with(tolower) # convert all column headers to lowercase # check values of Date column for mistakes...looks good unique(dat$date) # check values of score column for mistakes...one row to fix unique(dat$score) # one row contains no data for the score (row 767)...needs to be fixed dat %>% filter(score == "") # check values of rings_comb column for mistakes...need to remove the non-unique combos (rows with XB and XW) unique(dat$rings_comb) dat %>% filter(rings_comb %in% c("XB", "XW")) dat <- dat %>% filter(!grepl("XB", rings_comb) & !grepl("XW", rings_comb)) # import sex, migratory status, and age at banding information dat_sexes <- read.csv("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/sex_status_banding.csv", colClasses = "character") %>% mutate(Rings_comb = str_sub(Rings_comb, 1, 4)) %>% rename_with(tolower) %>% # convert all column headers to lowercase distinct() # checkout the sex-type data # number of birds with sex-type data dat_sexes %>% pull(rings_comb) %>% unique() %>% length() # identify duplicates...none, good dat_sexes[which(duplicated(dat_sexes)), ] # assess all birds with migratory status data dat_sexes %>% filter(migratory_status %in% c("R", "M")) # mutate the Date column into a date variable dat <- dat %>% mutate(date = paste(substring(dat$date, first = 7, last = 10), substring(date, first = 4, last = 5), substring(date, first = 1, last = 2), sep = "-") %>% as.Date()) %>% # subset to data with Molt == 1 filter(molt == 1) %>% # specify the season as the first calender year mutate(season = ifelse(month(date) < 7, year(date) - 1, year(date))) %>% # change to Julian date shifted for the Southern Hemisphere (1 = July 1) mutate(date_J = as.numeric(format(date + 181, "%j"))) %>% # join the sexes provided by Ailsa (note many-to-many is fine) left_join(., dat_sexes, by = "rings_comb", relationship = "many-to-many") %>% # remove individuals with unknown sex filter(sex != "U") %>% filter(sex != "N") %>% # add a ranking variable to sort the facets of the sampling distribution plots group_by(rings_comb) %>% mutate(n_photos = n(), n_scores = n_distinct(score)) %>% arrange(desc(n_scores)) %>% mutate( banding_date = paste(substring(banding_date, first = 7, last = 10), substring(banding_date, first = 4, last = 5), substring(banding_date, first = 1, last = 2), sep = "-") %>% as.Date(), season_ringed = ifelse(age_at_banding == "P" , ifelse(month(banding_date) < 7, year(banding_date) - 1, year(banding_date)), NA), age = ifelse(age_at_banding == "P", as.numeric(season - season_ringed), NA))# i should find a better way to calculate age # check max and min dates for each season dat %>% group_by(season) %>% summarise(n_obs = n(), min_date = min(date), max_date = max(date)) # summarize the number of individuals in the dataset...94 unique combos dat %>% pull(rings_comb) %>% unique() %>% length() # summarise the number of seasons for each individual...63 have more than one season of data dat %>% select(rings_comb, season) %>% distinct() %>% group_by(rings_comb) %>% summarise(n_seasons = n()) %>% arrange(desc(n_seasons)) %>% filter(n_seasons > 1) %>% nrow() # summarise the number of usable photos for each individual... dat %>% group_by(rings_comb, season) %>% summarise(n_photos = n()) %>% arrange(desc(n_photos)) ``` ### Plot the sampling distribution of Ailsa's photos across the seasons #### 2021-2022 season data ```{r, fig.height=9} # 2021-2022 season data ggplot(data = dat %>% filter(season == 2021 & sex == "F") %>% mutate(score = as.numeric(score), rings_comb_ = reorder(rings_comb, n_photos, decreasing = TRUE))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ rings_comb_, ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) + xlab("week") + ggtitle("Females from the 2021-2022 breeding season in Kaikoura") ``` ```{r, fig.height=8} ggplot(data = dat %>% filter(season == 2021 & sex == "M") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2021-07-05"), as.Date("2022-05-01"))) + xlab("week") + ggtitle("Males from the 2021-2022 breeding season in Kaikoura") ``` #### 2022-2023 season data ```{r, fig.height=10} # 2022-2023 season data (not as good coverage as the previous season) ggplot(data = dat %>% filter(season == 2022 & sex == "F") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) + xlab("week") + ggtitle("Females from the 2022-2023 breeding season in Kaikoura") ``` ```{r, fig.height=7} ggplot(data = dat %>% filter(season == 2022 & sex == "M") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2022-07-05"), as.Date("2023-05-01"))) + xlab("week") + ggtitle("Males from the 2022-2023 breeding season in Kaikoura") ``` #### 2023-2024 season data ```{r, fig.height=10} # 2023-2024 season data ggplot(data = dat %>% filter(season == 2023 & sex == "F") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2023-07-01"), as.Date("2024-05-01"))) + xlab("week") + ggtitle("Females from the 2023-2024 breeding season in Kaikoura") ``` ```{r, fig.height=6} ggplot(data = dat %>% filter(season == 2023 & sex == "M") %>% mutate(score = as.numeric(score))) + geom_point(aes(y = 1, x = date, fill = score), pch = 21, color = "black", size = 3) + facet_wrap(. ~ reorder(rings_comb, n_photos, decreasing = TRUE), ncol = 1, strip.position = "right") + scale_fill_gradient(high = "#cc4c02", low = "white") + theme_bw() + theme(axis.title.y = element_blank(), axis.text.y = element_blank(), panel.grid.major.y = element_blank(), panel.grid.minor = element_blank(), axis.ticks = element_blank(), strip.text.y.right = element_text(angle = 0)) + scale_x_date(date_labels = "%W", expand = c(0.01, 0.01), date_breaks = "3 week", limits = c(as.Date("2023-07-05"), as.Date("2024-05-01"))) + xlab("week") + ggtitle("Males from the 2023-2024 breeding season in Kaikoura") ``` ## Statistical analysis of moult dynamics ### Sex-specific differences in breeding plumage ```{r, include=FALSE} # extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual. ind_breeding_scores <- dat %>% mutate(score = as.numeric(score)) %>% mutate(breeding_season = ifelse(month(date) %in% c(8, 9, 10, 11, 12), 1, 0)) %>% filter(breeding_season == 1) %>% group_by(season, rings_comb, sex) %>% summarise(max_breeding_score = max(score)) %>% ungroup() # Assess sample sizes of each sex (49 females, 35 males) ind_breeding_scores %>% group_by(sex) %>% summarise(n_distinct(rings_comb)) # linear mixed model for the difference in max score between sexes mod_max_score <- lmer(max_breeding_score ~ sex + (1 | rings_comb) + (1 | season), data = ind_breeding_scores) # Derive confidence intervals of effect sizes from parametric bootstrapping tidy_mod_max_score <- tidy(mod_max_score, conf.int = TRUE, conf.method = "boot", nsim = 1000) tbl_regression(mod_max_score, intercept = TRUE, label = list(sex ~ "Sex")) # run rptR to obtain repeatabilities of random effects rpt_mod_max_score <- rpt(max_breeding_score ~ sex + (1 | rings_comb) + (1 | season), grname = c("rings_comb", "season", "Fixed"), data = ind_breeding_scores, datatype = "Gaussian", nboot = 1000, npermut = 1000, ratio = TRUE, adjusted = TRUE, ncores = 4, parallel = TRUE) # run partR2 on each model to obtain marginal R2, parameter estimates, and beta # weights R2m_mod_max_score <- partR2(mod_max_score, partvars = c("sex"), R2_type = "marginal", nboot = 1000, CI = 0.95, max_level = 1) R2c_mod_max_score <- partR2(mod_max_score, partvars = c("sex"), R2_type = "conditional", nboot = 1000, CI = 0.95, max_level = 1) stats_mod_max_score <- list(mod = mod_max_score, tidy = tidy_mod_max_score, rptR = rpt_mod_max_score, partR2m = R2m_mod_max_score, partR2c = R2c_mod_max_score, data = ind_breeding_scores) ``` ```{r, include=FALSE} # saveRDS(stats_mod_max_score, file = "out/stats_mod_max_score.rds") stats_mod_max_score <- readRDS(file = "/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/out/stats_mod_max_score.rds") ``` ```{r, include=FALSE} #### Table of effect sizes ---- # Retrieve sample sizes sample_sizes <- ind_breeding_scores %>% ungroup() %>% summarise(Year = n_distinct(season), Individual = n_distinct(rings_comb), Observations = nrow(.)) sample_sizes <- as.data.frame(t(as.data.frame(sample_sizes))) %>% rownames_to_column("term") %>% rename(estimate = V1) %>% mutate(stat = "n") # clean model component names mod_comp_names <- data.frame(comp_name = c("Male (Sex)", "Total Marginal \U1D479\U00B2", "Sex", "Total Conditional \U1D479\U00B2", "Individual", "Year", "Residual", "Individual", "Year", "Residual", "Years", "Individuals", "Observations")) # Fixed effect sizes (non-standardized) fixefTable <- stats_mod_max_score$tidy %>% dplyr::filter(effect == "fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "fixed") # Fixed effect sizes (standardized) fixef_bw_Table <- stats_mod_max_score$partR2m$BW %>% as.data.frame() %>% mutate(stat = "fixed_bw") %>% rename(conf.low = CI_lower, conf.high = CI_upper) # Semi-partial R2 estimates R2Table <- bind_rows(stats_mod_max_score$partR2m$R2, stats_mod_max_score$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>% as.data.frame() %>% mutate(stat = "partR2") %>% rename(conf.low = CI_lower, conf.high = CI_upper) # Random effects variances ranefTable <- stats_mod_max_score$tidy %>% dplyr::filter(effect == "ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "rand") %>% rename(term = group) %>% mutate(estimate = estimate^2, conf.high = conf.high^2, conf.low = conf.low^2) # Adjusted repeatabilities coefRptTable <- stats_mod_max_score$rptR$R_boot %>% dplyr::select(-Fixed) %>% mutate(residual = 1 - rowSums(.)) %>% apply(., 2, function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% t() %>% as.data.frame() %>% rownames_to_column("term") %>% rename(estimate = V1, conf.low = `2.5%`, conf.high = `97.5%`) %>% mutate(stat = "RptR") # Store all parameters into a single table and clean it up allCoefs_mod <- bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>% bind_cols(., mod_comp_names) %>% mutate(coefString = ifelse(!is.na(conf.low), paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"), NA), effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)), rep("Partitioned \U1D479\U00B2", nrow(R2Table)), rep("Random effects \U1D70E\U00B2", nrow(ranefTable)), rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)), rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything()) # draw gt table mod_max_score_table <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>% gt(rowname_col = "row", groupname_col = "effect") %>% cols_label(comp_name = html("Banded Dotterel rufous band score"), estimate = "Mean estimate", coefString = "95% confidence interval") %>% fmt_number(columns = c(estimate), rows = 1:10, decimals = 2, use_seps = FALSE) %>% fmt_number(columns = c(estimate), rows = 11:13, decimals = 0, use_seps = FALSE) %>% sub_missing(columns = 1:4, missing_text = "") %>% cols_align(align = "left", columns = c(comp_name)) %>% tab_options(row_group.font.weight = "bold", row_group.background.color = brewer.pal(9,"Greys")[3], table.font.size = 12, data_row.padding = 3, row_group.padding = 4, summary_row.padding = 2, column_labels.font.size = 14, row_group.font.size = 12, table.width = pct(80)) ``` #### sex-specific breeding plumage distributions ```{r} ind_breeding_scores %>% group_by(rings_comb, sex) %>% summarise(max_breeding_score = max(max_breeding_score)) %>% ungroup() %>% group_by(sex, max_breeding_score) %>% summarise(n = n()) %>% bind_rows(., data.frame(sex = c("F", "M"), max_breeding_score = c(7, 3), n = c(0, 0))) %>% ungroup() %>% ggplot() + geom_bar(aes(x = max_breeding_score, y = n, fill = sex), alpha = 0.5, stat = "identity", position = "dodge", width = 0.5) + luke_theme + theme(legend.position = c(0.15, 0.9)) + xlab("Maximum individual rufous band score during breeding season") + ylab("Frequency") + scale_colour_brewer(palette = "Dark2", direction = -1, name = "Sex", labels = c("Female (N = 49)", "Male (N = 35)")) + scale_fill_brewer(palette = "Dark2", direction = -1, name = "Sex", labels = c("Female (N = 49)", "Male (N = 35)")) ``` #### effect-size table for sex-specific breeding plumage Males tend to have a more full and immaculate rufous bands than females ```{r, include=FALSE} mod_max_score_table ``` ```{r, eval=FALSE} effects_mod_max_score <- as.data.frame(allEffects(mod_max_score)$sex) %>% mutate(x_mean = ifelse(sex == "F", as.numeric(factor(sex)) + 0.2, ifelse(sex == "M", as.numeric(factor(sex)) - 0.2, NA)), x_data = ifelse(sex == "F", as.numeric(factor(sex)), ifelse(sex == "M", as.numeric(factor(sex)), NA))) ggplot() + geom_errorbar(data = effects_mod_max_score, aes(x = x_mean, ymin = lower, ymax = upper), width = 0.05) + geom_point(data = effects_mod_max_score, aes(x = x_mean, y = fit, fill = sex), shape = 21, size = 5) + geom_jitter(data = ind_breeding_scores, aes(x = as.numeric(factor(sex)), y = max_breeding_score, fill = sex), width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75) + theme(legend.position = "none", text = element_text(family = "lato"), axis.ticks = element_blank(), axis.text = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), axis.text.y = element_text(size = 10, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), strip.background = element_rect(fill = 'grey90', color = NA), strip.text = element_text(size = 13, colour = "grey40"), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), axis.title.x = element_blank()) + scale_x_continuous(limits = c(0.8, 2.2), breaks = c(1.1, 1.9), labels = c("female", "male")) + ylab("maximum breeding plumage score observed") + scale_fill_brewer(palette = "Set2") + scale_y_continuous(limits = c(3.8, 7.2)) ``` ### Age- and sex-specific variation in breeding plumage (of known-aged birds) ```{r, include=FALSE} # extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual. ind_breeding_scores_age <- dat %>% mutate(score = as.numeric(score)) %>% mutate(breeding_season = ifelse(month(date) %in% c(8, 9, 10, 11, 12), 1, 0)) %>% filter(breeding_season == 1) %>% group_by(season, rings_comb, sex, age) %>% summarise(max_breeding_score = max(score)) %>% ungroup() %>% filter(!is.na(age)) # Assess sample sizes of each sex and season (14 females, 10 males) ind_breeding_scores_age %>% group_by(sex) %>% summarise(n_distinct(rings_comb)) # linear mixed model for the difference in max score between sexes across age mod_max_score_age_sex <- lmer(max_breeding_score ~ sex + age + (1 | rings_comb) + (1 | season), data = ind_breeding_scores_age) # Derive confidence intervals of effect sizes from parametric bootstrapping tidy_mod_max_score_age_sex <- tidy(mod_max_score_age_sex, conf.int = TRUE, conf.method = "boot", nsim = 1000) tbl_regression(mod_max_score_age_sex, intercept = TRUE, label = list(sex ~ "Sex", age ~ "Age")) # run rptR to obtain repeatabilities of random effects rpt_mod_max_score_age_sex <- rpt(max_breeding_score ~ sex + age + (1 | rings_comb) + (1 | season), grname = c("rings_comb", "season", "Fixed"), data = ind_breeding_scores_age, datatype = "Gaussian", nboot = 1000, npermut = 1000, ratio = TRUE, adjusted = TRUE, ncores = 4, parallel = TRUE) # run partR2 on each model to obtain marginal R2, parameter estimates, and beta # weights R2m_mod_max_score_age_sex <- partR2(mod_max_score_age_sex, partvars = c("sex", "age"), R2_type = "marginal", nboot = 1000, CI = 0.95, max_level = 1) R2c_mod_max_score_age_sex <- partR2(mod_max_score_age_sex, partvars = c("sex", "age"), R2_type = "conditional", nboot = 1000, CI = 0.95, max_level = 1) stats_mod_max_score_age_sex <- list(mod = mod_max_score_age_sex, tidy = tidy_mod_max_score_age_sex, rptR = rpt_mod_max_score_age_sex, partR2m = R2m_mod_max_score_age_sex, partR2c = R2c_mod_max_score_age_sex, data = ind_breeding_scores_age) ``` ```{r, include=FALSE} # saveRDS(stats_mod_max_score_age_sex, file = "out/stats_mod_max_score_age_sex.rds") stats_mod_max_score_age_sex <- readRDS(file = "/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/out/stats_mod_max_score_age_sex.rds") ``` ```{r, include=FALSE} #### Table of effect sizes ---- # Retrieve sample sizes sample_sizes <- ind_breeding_scores_age %>% ungroup() %>% summarise(Year = n_distinct(season), Individual = n_distinct(rings_comb), Observations = nrow(.)) sample_sizes <- as.data.frame(t(as.data.frame(sample_sizes))) %>% rownames_to_column("term") %>% rename(estimate = V1) %>% mutate(stat = "n") # clean model component names mod_comp_names <- data.frame(comp_name = c("Male (Sex)", "Age", "Total Marginal \U1D479\U00B2", "Sex", "Age", "Total Conditional \U1D479\U00B2", "Individual", "Year", "Residual", "Individual", "Year", "Residual", "Years", "Individuals", "Observations")) # Fixed effect sizes (non-standardized) fixefTable <- stats_mod_max_score_age_sex$tidy %>% dplyr::filter(effect == "fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "fixed") # Fixed effect sizes (standardized) fixef_bw_Table <- stats_mod_max_score_age_sex$partR2m$BW %>% as.data.frame() %>% mutate(stat = "fixed_bw") %>% rename(conf.low = CI_lower, conf.high = CI_upper) # Semi-partial R2 estimates R2Table <- bind_rows(stats_mod_max_score_age_sex$partR2m$R2, stats_mod_max_score_age_sex$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>% as.data.frame() %>% mutate(stat = "partR2") %>% rename(conf.low = CI_lower, conf.high = CI_upper) # Random effects variances ranefTable <- stats_mod_max_score_age_sex$tidy %>% dplyr::filter(effect == "ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "rand") %>% rename(term = group) %>% mutate(estimate = estimate^2, conf.high = conf.high^2, conf.low = conf.low^2) # Adjusted repeatabilities coefRptTable <- stats_mod_max_score_age_sex$rptR$R_boot %>% dplyr::select(-Fixed) %>% mutate(residual = 1 - rowSums(.)) %>% apply(., 2, function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% t() %>% as.data.frame() %>% rownames_to_column("term") %>% rename(estimate = V1, conf.low = `2.5%`, conf.high = `97.5%`) %>% mutate(stat = "RptR") # Store all parameters into a single table and clean it up allCoefs_mod <- bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>% bind_cols(., mod_comp_names) %>% mutate(coefString = ifelse(!is.na(conf.low), paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"), NA), effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)), rep("Partitioned \U1D479\U00B2", nrow(R2Table)), rep("Random effects \U1D70E\U00B2", nrow(ranefTable)), rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)), rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything()) # draw gt table mod_max_score_age_sex_table <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>% gt(rowname_col = "row", groupname_col = "effect") %>% cols_label(comp_name = html("Banded Dotterel rufous band score"), estimate = "Mean estimate", coefString = "95% confidence interval") %>% fmt_number(columns = c(estimate), rows = 1:12, decimals = 2, use_seps = FALSE) %>% fmt_number(columns = c(estimate), rows = 11:15, decimals = 0, use_seps = FALSE) %>% sub_missing(columns = 1:4, missing_text = "") %>% cols_align(align = "left", columns = c(comp_name)) %>% tab_options(row_group.font.weight = "bold", row_group.background.color = brewer.pal(9,"Greys")[3], table.font.size = 12, data_row.padding = 3, row_group.padding = 4, summary_row.padding = 2, column_labels.font.size = 14, row_group.font.size = 12, table.width = pct(80)) ``` #### effect-size table for sex- and age-specific variation in breeding plumage no effect of age on the maximum extent of an individual's rufous band score while controlling for underlying sex-specific variation ```{r, include=FALSE} mod_max_score_age_sex_table ``` ```{r, include=FALSE} # extract predicted trends pr <- ggeffects::predict_response(mod_max_score_age_sex) # join the back-transformed dates to model fits mod_max_score_age_sex_fits <- as.data.frame(pr) ggplot() + geom_ribbon(data = mod_max_score_age_sex_fits, aes(x = age.x, ymin = age.conf.low, ymax = age.conf.high), fill = "steelblue", alpha = 0.5) + geom_line(data = mod_max_score_age_sex_fits, aes(x = age.x, y = age.predicted), color = "steelblue") + geom_jitter(data = ind_breeding_scores_age, aes(x = age, y = max_breeding_score), width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75, fill = "steelblue", alpha = 1) + theme(legend.position = "none", text = element_text(family = "lato"), axis.ticks = element_blank(), axis.text = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), axis.text.y = element_text(size = 10, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), strip.background = element_rect(fill = 'grey90', color = NA), strip.text = element_text(size = 13, colour = "grey40"), panel.grid.major = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank()) + scale_x_continuous(limits = c(0.5, 6.5), breaks = c(1:6), labels = c(1:6)) + ylab("maximum breeding plumage score observed") + xlab("age") + scale_fill_brewer(palette = "Set2") + scale_y_continuous(limits = c(3.8, 7.2)) ``` ### Migratory status- and sex-specific variation in breeding plumage (of known-status birds) ```{r, include=FALSE} # extract the core breeding months (i.e., when presumably all birds are at their maximum breeding plumage), and determine the maximum score for each individual. ind_breeding_scores_status <- dat %>% mutate(score = as.numeric(score)) %>% mutate(breeding_season = ifelse(month(date) %in% c(8, 9, 10, 11, 12), 1, 0)) %>% filter(breeding_season == 1) %>% group_by(season, rings_comb, sex, age, migratory_status) %>% summarise(max_breeding_score = max(score)) %>% ungroup() %>% filter(migratory_status %in% c("M", "R")) # Assess sample sizes of each sex and season # (13 females (2 migrant, 11 resident), 8 males (3 migrant, 5 resident)) ind_breeding_scores_status %>% group_by(sex, migratory_status) %>% summarise(n_distinct(rings_comb)) # linear mixed model for the difference in max score between sexes mod_max_score_status_sex <- lmer(max_breeding_score ~ sex + migratory_status + (1 | rings_comb), data = ind_breeding_scores_status) # # detected singularity issue in a model with (1 | season), summary showed # that variance component of this random effect was zero. # # Model with one random effect mod1 <- lmer(max_breeding_score ~ sex + migratory_status + season + (1 | rings_comb), data = ind_breeding_scores_status) mod2 <- lmer(max_breeding_score ~ sex + migratory_status + (1 | season), data = ind_breeding_scores_status) mod3 <- lmer(max_breeding_score ~ sex + migratory_status + (1 | rings_comb), data = ind_breeding_scores_status) # Compare with the full model mod_full <- lmer(max_breeding_score ~ sex + migratory_status + (1 | rings_comb) + (1 | season), data = ind_breeding_scores_status) # Check the summaries summary(mod1) summary(mod2) summary(mod3) summary(mod_full) # Compare models using AIC AIC(mod1, mod2, mod3, mod_full) # conclude that the term (1 | season) is not needed for statistical reasons, but also since we are not interested in population-level annual variation in the max-breeding plumage scores we can drop this # Derive confidence intervals of effect sizes from parametric bootstrapping tidy_mod_max_score_status_sex <- tidy(mod_max_score_status_sex, conf.int = TRUE, conf.method = "boot", nsim = 1000) tbl_regression(mod_max_score_status_sex, intercept = TRUE, label = list(sex ~ "Sex", migratory_status ~ "Status")) # run rptR to obtain repeatabilities of random effects rpt_mod_max_score_status_sex <- rpt(max_breeding_score ~ sex + migratory_status + (1 | rings_comb), grname = c("rings_comb", "Fixed"), data = ind_breeding_scores_status, datatype = "Gaussian", nboot = 1000, npermut = 1000, ratio = TRUE, adjusted = TRUE, ncores = 4, parallel = TRUE) # run partR2 on each model to obtain marginal R2, parameter estimates, and beta # weights R2m_mod_max_score_status_sex <- partR2(mod_max_score_status_sex, partvars = c("sex", "migratory_status"), R2_type = "marginal", nboot = 1000, CI = 0.95, max_level = 1) R2c_mod_max_score_status_sex <- partR2(mod_max_score_status_sex, partvars = c("sex", "migratory_status"), R2_type = "conditional", nboot = 1000, CI = 0.95, max_level = 1) stats_mod_max_score_status_sex <- list(mod = mod_max_score_status_sex, tidy = tidy_mod_max_score_status_sex, rptR = rpt_mod_max_score_status_sex, partR2m = R2m_mod_max_score_status_sex, partR2c = R2c_mod_max_score_status_sex, data = ind_breeding_scores) ``` ```{r, include=FALSE} # saveRDS(stats_mod_max_score_status_sex, file = "out/stats_mod_max_score_status_sex.rds") stats_mod_max_score_status_sex <- readRDS(file = "/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/out/stats_mod_max_score_status_sex.rds") ``` ```{r, include=FALSE} #### Table of effect sizes ---- # Retrieve sample sizes sample_sizes <- ind_breeding_scores_status %>% ungroup() %>% summarise(Year = n_distinct(season), Individual = n_distinct(rings_comb), Observations = nrow(.)) sample_sizes <- as.data.frame(t(as.data.frame(sample_sizes))) %>% rownames_to_column("term") %>% rename(estimate = V1) %>% mutate(stat = "n") # clean model component names mod_comp_names <- data.frame(comp_name = c("Male (Sex)", "Resident (Migratory Status)", "Total Marginal \U1D479\U00B2", "Sex", "Migratory Status", "Total Conditional \U1D479\U00B2", "Individual", "Residual", "Individual", "Residual", "Years", "Individuals", "Observations")) # Fixed effect sizes (non-standardized) fixefTable <- stats_mod_max_score_status_sex$tidy %>% dplyr::filter(effect == "fixed") %>% dplyr::select(term, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "fixed") # Fixed effect sizes (standardized) fixef_bw_Table <- stats_mod_max_score_status_sex$partR2m$BW %>% as.data.frame() %>% mutate(stat = "fixed_bw") %>% rename(conf.low = CI_lower, conf.high = CI_upper) # Semi-partial R2 estimates R2Table <- bind_rows(stats_mod_max_score_status_sex$partR2m$R2, stats_mod_max_score_status_sex$partR2c$R2[1,]) %>% dplyr::select(term, estimate, CI_lower, CI_upper) %>% as.data.frame() %>% mutate(stat = "partR2") %>% rename(conf.low = CI_lower, conf.high = CI_upper) # Random effects variances ranefTable <- stats_mod_max_score_status_sex$tidy %>% dplyr::filter(effect == "ran_pars") %>% dplyr::select(group, estimate, conf.low, conf.high) %>% as.data.frame() %>% mutate(stat = "rand") %>% rename(term = group) %>% mutate(estimate = estimate^2, conf.high = conf.high^2, conf.low = conf.low^2) # Adjusted repeatabilities coefRptTable <- stats_mod_max_score_status_sex$rptR$R_boot %>% dplyr::select(-Fixed) %>% mutate(residual = 1 - rowSums(.)) %>% apply(., 2, function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% t() %>% as.data.frame() %>% rownames_to_column("term") %>% rename(estimate = V1, conf.low = `2.5%`, conf.high = `97.5%`) %>% mutate(stat = "RptR") # Store all parameters into a single table and clean it up allCoefs_mod <- bind_rows(fixef_bw_Table, R2Table, ranefTable, coefRptTable, sample_sizes) %>% bind_cols(., mod_comp_names) %>% mutate(coefString = ifelse(!is.na(conf.low), paste0("[", round(conf.low, 2), ", ", round(conf.high, 2), "]"), NA), effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)), rep("Partitioned \U1D479\U00B2", nrow(R2Table)), rep("Random effects \U1D70E\U00B2", nrow(ranefTable)), rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)), rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>% dplyr::select(effect, everything()) ``` ```{r} # draw gt table mod_max_score_table_status_sex <- allCoefs_mod %>% dplyr::select(effect, comp_name, estimate, coefString) %>% gt(rowname_col = "row", groupname_col = "effect") %>% cols_label(comp_name = html("Banded Dotterel rufous band score"), estimate = "Mean estimate", coefString = "95% confidence interval") %>% fmt_number(columns = c(estimate), rows = 1:10, decimals = 2, use_seps = FALSE) %>% fmt_number(columns = c(estimate), rows = 11:13, decimals = 0, use_seps = FALSE) %>% sub_missing(columns = 1:4, missing_text = "") %>% cols_align(align = "left", columns = c(comp_name)) %>% tab_options(row_group.font.weight = "bold", row_group.background.color = brewer.pal(9,"Greys")[3], table.font.size = 12, data_row.padding = 3, row_group.padding = 4, summary_row.padding = 2, column_labels.font.size = 14, row_group.font.size = 12, table.width = pct(80)) ``` #### effect-size table and plot for sex- and migratory-status breeding plumage no effect of migratory status on an individual's rufous band score while controlling for underlying sex-specific variation ```{r, include=FALSE} mod_max_score_table_status_sex effects_mod_max_score_status_sex2 <- as.data.frame(allEffects(mod_max_score_status_sex)$migratory_status) %>% mutate(x_mean = ifelse(migratory_status == "M", as.numeric(factor(migratory_status)) + 0.2, ifelse(migratory_status == "R", as.numeric(factor(migratory_status)) - 0.2, NA)), x_data = ifelse(migratory_status == "M", as.numeric(factor(migratory_status)), ifelse(migratory_status == "R", as.numeric(factor(migratory_status)), NA))) ggplot() + geom_errorbar(data = effects_mod_max_score_status_sex2, aes(x = x_mean, ymin = lower, ymax = upper), width = 0.05) + geom_point(data = effects_mod_max_score_status_sex2, aes(x = x_mean, y = fit, fill = migratory_status), shape = 21, size = 5) + geom_jitter(data = ind_breeding_scores_status, aes(x = as.numeric(factor(migratory_status)), y = max_breeding_score, fill = migratory_status), width = 0.05, height = 0.1, shape = 21, size = 3, alpha = 0.75) + theme(legend.position = "none", text = element_text(family = "lato"), axis.ticks = element_blank(), axis.text = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), axis.text.y = element_blank(), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), axis.title.x = element_blank()) + scale_x_continuous(limits = c(0.8, 2.2), breaks = c(1.1, 1.9), labels = c("migrant", "resident")) + ylab("maximum breeding plumage score observed") + scale_fill_brewer(palette = "Set1") + scale_y_continuous(limits = c(3.8, 7.2)) ``` ### Within-individual moult dynamics across sexes ```{r} # calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in sex while accounting for sex differences in maximum rufous band score. ind_prop_molt_scores <- dat %>% mutate(score = as.numeric(score)) %>% left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% filter(!is.na(max_breeding_score)) %>% mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>% select(rings_comb, season, date_J, sex, score, max_breeding_score, prop_molt_score) %>% arrange(rings_comb, season, date_J) %>% distinct() %>% filter(!is.na(prop_molt_score)) %>% group_by(sex) %>% mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>% ungroup() # Assess sample sizes of each sex ind_prop_molt_scores %>% group_by(sex) %>% summarise(n_distinct(rings_comb)) # mixed effects binomial model comparing sex and date effect on the changes in moult scores while controlling for between sex differences in the maximum breeding score mod1 <- lme4::glmer(prop_molt_score ~ date_J * sex + std_max_breeding_score + (1 | rings_comb), data = ind_prop_molt_scores, family = binomial, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))) # strong date effect, but no sex effect tbl_regression(mod1, intercept = TRUE, label = list(date_J ~ "Date", sex ~ "Sex")) # extract predicted trends pr <- ggeffects::predict_response(mod1, c("date_J [30:293]", "sex")) # back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365)) # join the back-transformed dates to model fits mod1_fits <- as.data.frame(pr) %>% rename(date_J = x, sex = group) %>% left_join(., dates_for_plot, by = "date_J") ind_prop_molt_scores <- ind_prop_molt_scores %>% left_join(., dates_for_plot, by = "date_J") # plot the model ggplot() + geom_line(data = mod1_fits, aes(x = date, y = predicted, color = sex)) + geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = sex), lwd = 1, alpha = 0.25) + geom_jitter(data = ind_prop_molt_scores, aes(x = date, y = prop_molt_score, color = sex), alpha = 0.25, shape = 20, width = 0, height = 0.025) + theme(legend.position = c(0.3, 0.2), legend.justification = c(1, 0), text = element_text(family = "lato"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + scale_colour_brewer(palette = "Dark2", direction = -1, name = "Sex", labels = c("Female (N = 49)", "Male (N = 35)")) + scale_fill_brewer(palette = "Dark2", direction = -1, name = "Sex", labels = c("Female (N = 49)", "Male (N = 35)")) ``` ```{r} # wrangle the dates of plumage retention and moult for subsequent analysis # Step 1: Calculate the max_score for each individual and season max_scores <- ind_prop_molt_scores %>% group_by(rings_comb, season) %>% summarise(max_score = max(prop_molt_score, na.rm = TRUE), .groups = "drop") # Step 2: Merge the max_scores back into the original data mid_point_dates <- ind_prop_molt_scores %>% left_join(max_scores, by = c("rings_comb", "season")) %>% # Add the max_score group_by(rings_comb, season) %>% arrange(date_J) %>% # Ensure data is sorted by date_J # Step 3: Find the latest date_J for the maximum prop_molt_score filter(prop_molt_score == max_score) %>% slice_max(order_by = date_J, n = 1, with_ties = FALSE) %>% rename(latest_max_date_J = date_J) %>% # Step 4: Find the earliest date_J where prop_molt_score < max_score right_join( ind_prop_molt_scores %>% left_join(max_scores, by = c("rings_comb", "season")) %>% # Add the max_score here too group_by(rings_comb, season) %>% arrange(date_J) %>% filter(prop_molt_score < max_score) %>% slice_min(order_by = date_J, n = 1, with_ties = FALSE) %>% rename(earliest_decline_date_J = date_J), by = c("rings_comb", "season") ) %>% # Step 5: Compute the mid-point date_J mutate(mid_point_date_J = round((latest_max_date_J + earliest_decline_date_J) / 2)) %>% # Select relevant columns select(rings_comb, season, latest_max_date_J, earliest_decline_date_J, mid_point_date_J) %>% distinct() dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365)) moult_commence_ring_season <- left_join(mid_point_dates, dates_for_plot %>% rename(mid_point_date_J = date_J), by = "mid_point_date_J") ``` ### Within-individual moult dynamics across migratory status ```{r} # calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation in migratory status. ind_prop_molt_scores <- dat %>% mutate(score = as.numeric(score)) %>% left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% filter(!is.na(max_breeding_score)) %>% mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>% select(rings_comb, season, date_J, sex, migratory_status, score, max_breeding_score, prop_molt_score) %>% arrange(rings_comb, season, date_J) %>% distinct() %>% filter(!is.na(prop_molt_score)) %>% filter(migratory_status %in% c("R", "M")) # Assess sample sizes of each migratory status ind_prop_molt_scores %>% group_by(migratory_status) %>% summarise(n_distinct(rings_comb)) # mixed effects binomial model comparing migratory status and date effect on the changes in moult scores mod1_status <- lme4::glmer(prop_molt_score ~ date_J * migratory_status + (1 | rings_comb), data = ind_prop_molt_scores, family = binomial, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))) # strong date effect, but no status effect tbl_regression(mod1_status, intercept = TRUE, label = list(date_J ~ "Date", migratory_status ~ "Migratory status")) # extract predicted trends pr <- ggeffects::predict_response(mod1_status, c("date_J [30:293]", "migratory_status")) # back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365)) # join the back-transformed dates to model fits mod1_fits <- as.data.frame(pr) %>% rename(date_J = x, migratory_status = group) %>% left_join(., dates_for_plot, by = "date_J") ind_prop_molt_scores <- ind_prop_molt_scores %>% left_join(., dates_for_plot, by = "date_J") # plot the model ggplot() + geom_line(data = mod1_fits, aes(x = date, y = predicted, color = migratory_status)) + geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = migratory_status), lwd = 1, alpha = 0.25) + geom_jitter(data = ind_prop_molt_scores, aes(x = date, y = prop_molt_score, color = migratory_status), alpha = 0.25, shape = 20, width = 0, height = 0.025) + theme(legend.position = c(0.3, 0.2), legend.justification = c(1, 0), text = element_text(family = "lato"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month") + scale_colour_brewer(palette = "Set1", direction = -1, name = "Migratory status", labels = c("Resident (N = 16)", "Migratory (N = 5)")) + scale_fill_brewer(palette = "Set1", direction = -1, name = "Migratory status", labels = c("Resident (N = 16)", "Migratory (N = 5)")) ``` ### Within-individual moult dynamics across ages ```{r} # calculate the individual proportional moult scores by comparing each score to a given individual's max (determined in previous chunk). Piece apart variation between age groups. ind_prop_molt_scores <- dat %>% mutate(score = as.numeric(score)) %>% left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% filter(!is.na(max_breeding_score)) %>% mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>% select(rings_comb, season, date_J, sex, age, score, max_breeding_score, prop_molt_score) %>% arrange(rings_comb, season, date_J) %>% distinct() %>% filter(!is.na(prop_molt_score)) %>% filter(!is.na(age)) # Assess sample sizes of each sex ind_prop_molt_scores %>% group_by(age) %>% summarise(n_distinct(rings_comb)) # mixed effects binomial model comparing sex and date effect on the changes in moult scores mod1_age <- lme4::glmer(prop_molt_score ~ date_J + age + # date_J * sex + (1 | rings_comb),# + (1 | season), data = ind_prop_molt_scores, family = binomial)#, # control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+06))) plot(allEffects(mod1_age)) # strong date effect, but no status effect tbl_regression(mod1_age, intercept = TRUE, label = list(date_J ~ "Date", age ~ "Age")) # extract predicted trends pr <- ggeffects::predict_response(mod1_age, c("date_J [30:293]", "age")) # back-transform the dates (i.e., Julian dates were optimized for Austral summer), note that the year is irrelevent here and is only kept for simplicity. We are mainly interested in matching the month and day part of the date string with the transformed Julian date dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), date_J = c(1:365)) # join the back-transformed dates to model fits mod1_fits <- as.data.frame(pr) %>% rename(date_J = x, age = group) %>% left_join(., dates_for_plot, by = "date_J") ind_prop_molt_scores <- ind_prop_molt_scores %>% left_join(., dates_for_plot, by = "date_J") # plot the model ggplot() + geom_line(data = mod1_fits, aes(x = date, y = predicted, color = age)) + geom_ribbon(data = mod1_fits, aes(x = date, ymax = conf.high, ymin = conf.low, fill = age), lwd = 1, alpha = 0.25) + geom_jitter(data = ind_prop_molt_scores, aes(x = date, y = prop_molt_score, color = as.factor(age)), alpha = 0.25, shape = 20, width = 0, height = 0.025) + theme(legend.position = c(0.3, 0.2), legend.justification = c(1, 0), text = element_text(family = "lato"), axis.ticks = element_blank(), axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1), axis.text.y = element_text(size = 10, colour = "grey40"), axis.title.y = element_text(size = 11, colour = "grey40"), plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"), panel.spacing = unit(1, "cm"), panel.background = element_rect(fill = 'transparent', color = NA), plot.background = element_rect(fill = 'transparent', color = NA), panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"), panel.grid.minor = element_blank(), panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"), axis.title.x = element_blank()) + ylab("proportion intact of individual-based maximum rufous band") + scale_x_date(date_labels = "%B", expand = c(0.01, 0.01), date_breaks = "1 month", limits = c(as.Date("2021-10-01"), as.Date("2022-05-01"))) + scale_colour_brewer(palette = "Spectral", direction = -1, name = "Age", labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) + scale_fill_brewer(palette = "Spectral", direction = -1, name = "Age", labels = c("1 (N = 6)", "2 (N = 5)", "3 (N = 12)", "4 (N = 13", "5 (N = 4)", "6 (N = 2)")) ``` ## Breeding status analysis ```{r} # explore datasets provided by Ted ## 2021 season # import and consolidate key columns breeding_data_2021 <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx", sheet = "Kaikoura_NestVisit2021_1", col_types = "text") %>% mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% mutate(date = as.Date(visit_date_nz)) %>% select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% rename(nest_fate = `Nest Status`, eggs = `Egg count`, hatched = `Number Hatched`, fledged = `Number Fledged`, chicks = `Number of chicks seen`) # Define a function to extract the desired patterns and ensure enough columns extract_patterns <- function(text) { # Replace vertical bar '|' with an empty string text <- gsub("\\|", "", text) # Extract all 4 or 5-character long texts that start with R or r matches_r <- str_extract_all(text, "\\b[Rr]\\w{3,4}\\b")[[1]] # Extract UB, UN, UBF, UBM as standalone patterns matches_ub_un <- str_extract_all(text, "\\b[Uu][BbNnFfMm]\\b")[[1]] # Combine matches matches <- c(matches_r, matches_ub_un) # Ensure all matches are capitalized matches <- toupper(matches) # Ensure there are exactly two columns, filling with NA if necessary result <- c(matches, rep(NA, 2 - length(matches))) # Return the result as a named vector return(setNames(result, c("parent1", "parent2"))) } # Apply the function to the text column and store the results in new columns extracted_parents <- t(sapply(breeding_data_2021$NestID, extract_patterns)) breeding_data_2021 <- bind_cols(breeding_data_2021, extracted_parents) %>% select(-parent2) %>% rename(parent = parent1) %>% bind_rows(bind_cols(breeding_data_2021, extracted_parents) %>% select(-parent1) %>% rename(parent = parent2)) %>% filter(!is.na(parent)) %>% filter(parent %in% (dat %>% filter(season == 2021) %>% pull(rings_comb) %>% unique())) %>% filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>% group_by(parent) %>% arrange(desc(date)) %>% slice(1) %>% ungroup() %>% arrange(desc(date)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>% rename(last_date_breeding = date) %>% select(parent, last_date_breeding, nest_id, OBJECTID) # import and consolidate key columns breeding_data_2021_ <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_2021_Feb2022_upd_Nov22_Locations.xlsx", sheet = "Kaikoura_DotterelNest2021_0", col_types = "text") %>% mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"), fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"), found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"), fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"), found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% mutate(hatch_date = as.Date(hatch_date_nz), fail_date = as.Date(fail_date_nz), found_date = as.Date(found_date_nz)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% select(OBJECTID, nest_id, found_date, hatch_date, fail_date) breeding_2021_final <- left_join(breeding_data_2021, breeding_data_2021_, by = "nest_id") %>% mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | ((!is.na(hatch_date) & found_date < hatch_date)) | ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) | ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% group_by(parent) %>% mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% arrange(desc(days_diff)) %>% mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>% arrange(desc(last_date_breeding_final)) %>% select(parent, last_date_breeding_final) %>% # specify the season as the first calender year mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final))) ########### ## 2022 season # import and consolidate key columns breeding_data_2022 <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx", sheet = "Kaikoura_NestVisit_1", col_types = "text") %>% mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% mutate(date = as.Date(visit_date_nz)) %>% select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% rename(nest_fate = `Nest Status`, eggs = `Egg count`, hatched = `Number Hatched`, fledged = `Number Fledged`, chicks = `Number of chicks seen`) # Apply the function to the text column and store the results in new columns extracted_parents2 <- t(sapply(breeding_data_2022$NestID, extract_patterns)) breeding_data_2022 <- bind_cols(breeding_data_2022, extracted_parents2) %>% select(-parent2) %>% rename(parent = parent1) %>% bind_rows(bind_cols(breeding_data_2022, extracted_parents2) %>% select(-parent1) %>% rename(parent = parent2)) %>% filter(!is.na(parent)) %>% filter(parent %in% (dat %>% filter(season == 2022) %>% pull(rings_comb) %>% unique())) %>% filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>% group_by(parent) %>% arrange(desc(date)) %>% slice(1) %>% ungroup() %>% arrange(desc(date)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>% rename(last_date_breeding = date) %>% select(parent, last_date_breeding, nest_id, OBJECTID) # import and consolidate key columns breeding_data_2022_ <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Jan23_upd.xlsx", sheet = "Kaikoura_DotterelNest_0", col_types = "text") %>% mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"), fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"), found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"), fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"), found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% mutate(hatch_date = as.Date(hatch_date_nz), fail_date = as.Date(fail_date_nz), found_date = as.Date(found_date_nz)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% select(OBJECTID, nest_id, found_date, hatch_date, fail_date) breeding_2022_final <- left_join(breeding_data_2022, breeding_data_2022_, by = "nest_id") %>% mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | ((!is.na(hatch_date) & found_date < hatch_date)) | ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) | ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% group_by(parent) %>% mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% arrange(desc(days_diff)) %>% mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>% arrange(desc(last_date_breeding_final)) %>% select(parent, last_date_breeding_final) %>% # specify the season as the first calender year mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final))) ########### ## 2023 season # import and consolidate key columns breeding_data_2023 <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx", sheet = "Kaikoura_NestVisit_1", col_types = "text") %>% mutate(visit_date_ = as.POSIXct(as.numeric(`Visit Date`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(visit_date_nz = with_tz(visit_date_, tzone = "Pacific/Auckland")) %>% mutate(date = as.Date(visit_date_nz)) %>% select(OBJECTID, NestID, date, `Nest Status`, `Egg count`, `Number Hatched`, Notes, `Number Fledged`, `Number of chicks seen`) %>% rename(nest_fate = `Nest Status`, eggs = `Egg count`, hatched = `Number Hatched`, fledged = `Number Fledged`, chicks = `Number of chicks seen`) # Apply the function to the text column and store the results in new columns extracted_parents3 <- t(sapply(breeding_data_2023$NestID, extract_patterns)) breeding_data_2023 <- bind_cols(breeding_data_2023, extracted_parents3) %>% select(-parent2) %>% rename(parent = parent1) %>% bind_rows(bind_cols(breeding_data_2023, extracted_parents3) %>% select(-parent1) %>% rename(parent = parent2)) %>% filter(!is.na(parent)) %>% filter(parent %in% (dat %>% filter(season == 2023) %>% pull(rings_comb) %>% unique())) %>% filter(nest_fate == "Occupied" | as.numeric(chicks) > 0) %>% group_by(parent) %>% arrange(desc(date)) %>% slice(1) %>% ungroup() %>% arrange(desc(date)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", NestID)) %>% rename(last_date_breeding = date) %>% select(parent, last_date_breeding, nest_id, OBJECTID) # import and consolidate key columns breeding_data_2023_ <- read_excel("/Users/luketheduke2/ownCloud/kemp_projects/bdot/moult/data/Kaikoura_Dotterel_Feb2024_upd.xlsx", sheet = "Kaikoura_DotterelNest_0", col_types = "text") %>% mutate(hatch_date_ = as.POSIXct(as.numeric(`Date Hatched`) * 86400, origin = "1899-12-30", tz = "UTC"), fail_date_ = as.POSIXct(as.numeric(`Date Failed`) * 86400, origin = "1899-12-30", tz = "UTC"), found_date_ = as.POSIXct(as.numeric(`Date Found`) * 86400, origin = "1899-12-30", tz = "UTC")) %>% mutate(hatch_date_nz = with_tz(hatch_date_, tzone = "Pacific/Auckland"), fail_date_nz = with_tz(fail_date_, tzone = "Pacific/Auckland"), found_date_nz = with_tz(found_date_, tzone = "Pacific/Auckland")) %>% mutate(hatch_date = as.Date(hatch_date_nz), fail_date = as.Date(fail_date_nz), found_date = as.Date(found_date_nz)) %>% mutate(nest_id = sub("^([^[:space:]]+).*", "\\1", `Nest ID`)) %>% select(OBJECTID, nest_id, found_date, hatch_date, fail_date) breeding_2023_final <- left_join(breeding_data_2023, breeding_data_2023_, by = "nest_id") %>% mutate(date_check = ifelse((!is.na(found_date) & is.na(fail_date) & is.na(hatch_date)) | ((!is.na(hatch_date) & found_date < hatch_date)) | ((!is.na(fail_date) & !is.na(hatch_date)) & hatch_date < fail_date) | ((!is.na(fail_date) & found_date < fail_date)), 1, 0)) %>% group_by(parent) %>% mutate(last_date_breeding_ = max(found_date, hatch_date, fail_date, na.rm = TRUE)) %>% mutate(days_diff = last_date_breeding_ - last_date_breeding) %>% arrange(desc(days_diff)) %>% mutate(last_date_breeding_final = as.Date(ifelse(days_diff > 0, last_date_breeding + ceiling(days_diff/2), last_date_breeding))) %>% arrange(desc(last_date_breeding_final)) %>% select(parent, last_date_breeding_final) %>% # specify the season as the first calender year mutate(season = ifelse(month(last_date_breeding_final) < 7, year(last_date_breeding_final) - 1, year(last_date_breeding_final))) # combine breeding data from the 3 seasons into one table breeding_data_all <- bind_rows(breeding_2021_final,breeding_2022_final,breeding_2023_final) %>% rename(rings_comb = parent) %>% # change to Julian date shifted for the Southern Hemisphere (1 = July 1) mutate(breeding_date_J = as.numeric(format(last_date_breeding_final + 181, "%j"))) %>% ungroup() %>% # standardize the last breeding date by year mutate(breeding_date_J_s = scale_by(breeding_date_J ~ season, ., scale = 0)) %>% # make the scaled date variable numeric class mutate(breeding_date_J_snum = as.numeric(breeding_date_J_s)) ``` ### Relationships between breeding schedule and plumage ```{r} ind_prop_molt_scores_breeding_date <- dat %>% mutate(score = as.numeric(score)) %>% left_join(., select(ind_breeding_scores, -sex), by = c("rings_comb","season")) %>% filter(!is.na(max_breeding_score)) %>% mutate(prop_molt_score = (score-1)/(max_breeding_score-1)) %>% select(rings_comb, season, date_J, sex, migratory_status, age, score, max_breeding_score, prop_molt_score) %>% arrange(rings_comb, season, date_J) %>% distinct() %>% filter(!is.na(prop_molt_score)) %>% full_join(breeding_data_all, by = c("rings_comb", "season")) %>% mutate(date_diff = abs(breeding_date_J - date_J)) %>% arrange(rings_comb, season, date_diff) %>% group_by(rings_comb, season) %>% slice(1) %>% filter(date_diff < 30) %>% dplyr::select(rings_comb, season, sex, breeding_date_J, breeding_date_J_snum, prop_molt_score, max_breeding_score, age) %>% distinct() %>% group_by(sex) %>% mutate(std_max_breeding_score = (max_breeding_score - mean(max_breeding_score, na.rm = TRUE)) / sd(max_breeding_score, na.rm = TRUE)) %>% ungroup() %>% mutate(breeding_date_J_s = scale(breeding_date_J)) %>% left_join(., moult_commence_ring_season, by = c("rings_comb", "season")) ``` #### effect of breeding plumage score on prolonged breeding do individuals with a fuller breeding plumage (for their sex) breed later than the rest of the population? No effect ```{r} mod_max_score_breeding_date <- lme4::lmer(breeding_date_J_snum ~ std_max_breeding_score + (1 | rings_comb), data = ind_prop_molt_scores_breeding_date) tbl_regression(mod_max_score_breeding_date, intercept = TRUE) plot(allEffects(mod_max_score_breeding_date), type = "response") ``` #### effect of breeding plumage score on prolonged breeding do individuals with a late breeding attempts tend to have a fuller breeding plumage? No effect ```{r} mod_breeding_date_max_score <- lme4::lmer(std_max_breeding_score ~ breeding_date_J_snum + (1 | rings_comb), data = ind_prop_molt_scores_breeding_date) summary(mod_max_score_breeding_date) plot(allEffects(mod_breeding_date_max_score), type = "response") ``` ### Relationships between breeding schedule and moult #### effect of prolonged breeding on moult do individuals with late breeding attempts have a delayed moult? Yes: individuals that breed later than the rest of the population in a given year tend to commence their complete/post-breeding moult later ```{r} mod_breeding_date_prop_molt_mid <- lme4::lmer(mid_point_date_J ~ breeding_date_J_snum + (1 | rings_comb), data = ind_prop_molt_scores_breeding_date) tbl_regression(mod_breeding_date_prop_molt_mid, intercept = TRUE) plot(allEffects(mod_breeding_date_prop_molt_mid), type = "response") ``` #### effect of prolonged breeding on prolonged retention of full plumage do individuals with late breeding attempts retain their full breeding plumage longer? Yes: individuals that breed later than the rest of the population in a given year tend to retain their full breeding plumage longer ```{r} mod_breeding_date_prop_molt_latest <- lm(latest_max_date_J ~ breeding_date_J_snum, data = ind_prop_molt_scores_breeding_date) tbl_regression(mod_breeding_date_prop_molt_latest, intercept = TRUE) plot(allEffects(mod_breeding_date_prop_molt_latest), type = "response") ``` ### Relationships between age and extended breeding period ```{r} # wrangle data to get ages and dates of last breeding for each season of all # birds available age_breeding_date <- dat %>% select(rings_comb, season, age, sex) %>% full_join(breeding_data_all, by = c("rings_comb", "season")) %>% distinct() ``` #### effect of age on prolonged breeding season do older individuals conclude their breeding season later? No: there is strong evidence that younger individuals tend to have the latest breeding attempts in the population for a given year. ```{r} mod_age_breeding_date <- lme4::lmer(breeding_date_J_snum ~ age + (1 | rings_comb), data = age_breeding_date) summary(mod_age_breeding_date) plot(allEffects(mod_age_breeding_date), type = "response") # strong date effect, but no status effect tbl_regression(mod_age_breeding_date, intercept = TRUE, label = list(age ~ "Age")) # extract predicted trends # I tried the whole range for the breeding timing (breeding_date_J [72:221]), but then decided to take some values with 20 days interval to visualize it better pr <- ggeffects::predict_response(mod_age_breeding_date, c("age")) dates_for_plot <- data.frame(date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), breeding_date_ = c(1:365)) CL_dates_for_plot <- data.frame(conf.low_date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), conf.low_ = c(1:365)) CH_dates_for_plot <- data.frame(conf.high_date = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), conf.high_ = c(1:365)) # join the back-transformed dates to model fits mod1_fits <- as.data.frame(pr) %>% rename(age = x, # sex = group, breeding_date = predicted) %>% mutate(breeding_date_ = round(breeding_date), conf.low_ = round(conf.low), conf.high_ = round(conf.high)) %>% left_join(., dates_for_plot, by = "breeding_date_") %>% left_join(., CL_dates_for_plot, by = "conf.low_") %>% left_join(., CH_dates_for_plot, by = "conf.high_") # plot the model ggplot() + geom_line(data = mod1_fits, aes(x = age, y = date)) + geom_ribbon(data = mod1_fits, aes(x = age, ymax = conf.high_date, ymin = conf.low_date), lwd = 1, alpha = 0.25) + luke_theme + # theme(legend.position = c(0.3, 0.2), # legend.justification = c(1, 0), # strip.background = element_blank(), # axis.title.x = element_blank(), # axis.text.x = element_text(size = 10, # angle = 45, # hjust = 1, # vjust = 1)) + ylab("date of last breeding") ``` ## Meeting notes: ### 13-Nov-24: Bart, Luke, and Bashar via zoom - add the breeding data to the analysis: for each individual include the date of last breeding (i.e., last date seen with nest or brood) as a predictor of molt timing. Also do an analysis relating the timing of last breeding to the age of the individual. Bashar will send Luke the breeding data he received from Ted. - drop pre-October data from the molt timing analysis. - next steps: Bashar to eventually write up this study as a report (which could eventually be used as his MSc thesis if he wants). Priority will be on him writing up a proposal for his thesis (which has a concrete deadline and is graded). Bashar will present the study at our weekly Monday meeting in the first few months of 2025.